# from python_speech_features import mfcc
# from python_speech_features import delta
# import numpy as np

# def get_feature(fs, signal):
#     remove_ixs = np.where(signal <= 0.0001)[0]
#     signal = np.delete(signal, remove_ixs)
#     mfcc_feature = mfcc(signal, fs)#signal - 需要用来计算特征的音频信号，应该是一个N*1的数组
#     if len(mfcc_feature) == 0:
#         print >> sys.stderr, "ERROR.. failed to extract mfcc feature:", len(signal)
#     mfcc_feature = delta(mfcc_feature,2)
#     return mfcc_feature
from __future__ import division
import numpy as np
from python_speech_features import sigproc
from scipy.fftpack import dct
import math
import wave


# def read(data_path):
#     '''读取语音信号
#     '''
#     # feat=[]
#     # appendEnergy = True
#     # if appendEnergy:
#     #     feat[:, 0] = np.log(energy)
#     wavepath = data_path
#     f = wave.open(wavepath, 'rb')
#     params = f.getparams()
#     nchannels, sampwidth, framerate, nframes = params[:4]  # 声道数、量化位数、采样频率、采样点数
#     str_data = f.readframes(nframes)  # 读取音频，字符串格式
#     f.close()
#     wavedata = np.fromstring(str_data, dtype=np.short)  # 将字符串转化为浮点型数据
#     wavedata = wavedata * 1.0 / (max(abs(wavedata)))  # wave幅值归一化
#     return wavedata, nframes, framerate


def enframe(data, win, inc):
    '''对语音数据进行分帧处理
    input:data(一维array):语音信号
      wlen(int):滑动窗长
      inc(int):窗口每次移动的长度
    output:f(二维array)每次滑动窗内的数据组成的二维array
    '''
    nx = len(data)  # 语音信号的长度
    try:
        nwin = len(win)
    except Exception as err:
        nwin = 1
    if nwin == 1:
        wlen = win
    else:
        wlen = nwin
    nf = int(np.fix((nx - wlen) / inc) + 1)  # 窗口移动的次数
    f = np.zeros((nf, wlen))  # 初始化二维数组
    indf = [inc * j for j in range(nf)]
    indf = (np.mat(indf)).T
    inds = np.mat(range(wlen))
    indf_tile = np.tile(indf, wlen)
    inds_tile = np.tile(inds, (nf, 1))
    mix_tile = indf_tile + inds_tile
    f = np.zeros((nf, wlen))
    for i in range(nf):
        for j in range(wlen):
            f[i, j] = data[mix_tile[i, j]]
    #print("fenzhen", len(f))
    #print("fenzhen", len(f[0]))
    return f


def point_check(wavedata, win, inc):
    '''语音信号端点检测
    input:wavedata(一维array)：原始语音信号
    output:StartPoint(int):起始端点
      EndPoint(int):终止端点
    '''
    # 1.计算短时过零率
    FrameTemp1 = enframe(wavedata[0:-1], win, inc)
    #frames = sigproc.framesig(signal, winlen*samplerate, winstep*samplerate, winfunc)#返回帧数组467*400分帧
    FrameTemp2 = enframe(wavedata[1:], win, inc)
    signs = np.sign(np.multiply(FrameTemp1, FrameTemp2))  # 计算每一位与其相邻的数据是否异号，异号则过零
    #sign()将元素值  >0 为1，小于0 为-1，等于0 为0，[-2.3,0,2.5]》》[-1,0,1]
    #np.multiply对应的是元素相乘a = np.mat([[1, 2, 3,]])，b = np.mat([[1, 2, 3,]])
    #c = np.multiply(a,b) = 1 4  9
    #将帧内所有sample平移1，再对应点做乘积，
    # 符号为负的则说明此处过零，只需将帧内所有负数乘积数目求出则得到该帧的过零率。
    signs = list(map(lambda x: [[i, 0][i > 0] for i in x], signs))
    #print(signs)
    signs = list(map(lambda x: [[i, 1][i < 0] for i in x], signs))
    #print("<0",signs)
    diffs = np.sign(abs(FrameTemp1 - FrameTemp2) - 0.01)
    diffs = list(map(lambda x: [[i, 0][i < 0] for i in x], diffs))
    zcr = list((np.multiply(signs, diffs)).sum(axis=1))
    #print("过零率",zcr)

    # 2.计算短时能量
    amp = list((abs(enframe(wavedata, win, inc))).sum(axis=1))
    # # 设置门限
    #print('设置di门限',round(np.mean(zcr) * 0.1))#np.mean()求取均值
    ZcrLow = max([round(np.mean(zcr) * 0.1), 3])  # 过零率低门限
    #print(ZcrLow)
    ZcrHigh = max([round(max(zcr) * 0.1), 5])  # 过零率高门限
    AmpLow = min([min(amp) * 10, np.mean(amp) * 0.2, max(amp) * 0.1])  # 能量低门限
    AmpHigh = max([min(amp) * 10, np.mean(amp) * 0.2, max(amp) * 0.1])  # 能量高门限
    # 端点检测
    MaxSilence = 8  # 最长语音间隙时间
    MinAudio = 16  # 最短语音时间
    Status = 0  # 状态0:静音段,1:过渡段,2:语音段,3:结束段
    HoldTime = 0  # 语音持续时间
    SilenceTime = 0  # 语音间隙时间
    #print('开始端点检测')
    StartPoint = 0
    for n in range(len(zcr)):
        if Status == 0 or Status == 1:
            if amp[n] > AmpHigh or zcr[n] > ZcrHigh:
                StartPoint = n - HoldTime
                Status = 2
                HoldTime = HoldTime + 1
                SilenceTime = 0
            elif amp[n] > AmpLow or zcr[n] > ZcrLow:
                Status = 1
                HoldTime = HoldTime + 1
            else:
                Status = 0
                HoldTime = 0
        elif Status == 2:
            if amp[n] > AmpLow or zcr[n] > ZcrLow:
                HoldTime = HoldTime + 1
            else:
                SilenceTime = SilenceTime + 1
                if SilenceTime < MaxSilence:
                    HoldTime = HoldTime + 1
                elif (HoldTime - SilenceTime) < MinAudio:
                    Status = 0
                    HoldTime = 0
                    SilenceTime = 0
                else:
                    Status = 3
        elif Status == 3:
            break
        if Status == 3:
            break
    HoldTime = HoldTime - SilenceTime
    EndPoint = StartPoint + HoldTime
    return FrameTemp1[StartPoint:EndPoint]


def get_feature(samplerate,signal,winlen=0.025,winstep=0.01,numcep=13,
         nfilt=26,nfft=512,lowfreq=0,highfreq=None,preemph=0.97,ceplifter=22,appendEnergy=True,
         winfunc=lambda x:np.ones((x,))):
    feat,energy = fbank(signal,samplerate,winlen,winstep,nfilt,nfft,lowfreq,highfreq,preemph,winfunc)
    feat = np.log(feat)
    feat = dct(feat, type=2, axis=1, norm='ortho')[:,:numcep]#[:,:numcep]','前面代表行，后面代表列
                                                              #numcep=13时候，代表所有行前13个元素
    feat = lifter(feat,ceplifter)
    if appendEnergy: feat[:,0] = np.log(energy) # replace first cepstral coefficient with log of frame energy
    feat = delta(feat,2)
    return feat
def fbank(signal,samplerate=16000,winlen=0.025,winstep=0.01,
          nfilt=26,nfft=512,lowfreq=0,highfreq=None,preemph=0.97,
          winfunc=lambda x:np.ones((x,))):
    """Compute Mel-filterbank energy features from an audio signal.

    :param signal: the audio signal from which to compute features. Should be an N*1 array
    :param samplerate: the samplerate of the signal we are working with.
    :param winlen: the length of the analysis window in seconds. Default is 0.025s (25 milliseconds)
    :param winstep: the step between successive windows in seconds. Default is 0.01s (10 milliseconds)
    :param nfilt: the number of filters in the filterbank, default 26.
    :param nfft: the FFT size. Default is 512.
    :param lowfreq: lowest band edge of mel filters. In Hz, default is 0.
    :param highfreq: highest band edge of mel filters. In Hz, default is samplerate/2
    :param preemph: apply preemphasis filter with preemph as coefficient. 0 is no filter. Default is 0.97.
    :param winfunc: the analysis window to apply to each frame. By default no window is applied. You can use np window functions here e.g. winfunc=np.hamming
    :returns: 2 values. The first is a np array of size (NUMFRAMES by nfilt) containing features. Each row holds 1 feature vector. The
        second return value is the energy in each frame (total energy, unwindowed)
    """
    highfreq= highfreq or samplerate/2#如果highfreq有输入则等于输入值，否则默认samplerate/2
    # print("语音信号",signal[-1]-signal[-1]*0.95)
    # print("语音信号", signal[-1])
    signal = sigproc.preemphasis(signal,preemph)#进行预加重处理，预加重系数为0.95将读取的语音信号输出滤波信号
    #print("滤波信号",signal[-1])
    #frames = sigproc.framesig(signal, winlen*samplerate, winstep*samplerate, winfunc)#返回帧数组467*400分帧
    frames = point_check(signal, nfft, inc=80)
    pspec = sigproc.powspec(frames,nfft)#计算每一帧的功率谱，如果帧是nxd矩阵，则输出为nx(nfft/2 1)。467*(256+1)
    energy = np.sum(pspec,1) # this stores the total energy in each frame这个装置储存了每一帧的总能量
    #np(a,axis)a:向量，数组或者矩阵，axis：可为None、整数、元祖，为0 时，压缩行，每一列元素相加，压缩成一行
    #print("能量",len(energy))                                                  #为1时，压缩列，每一行元素相加，压缩成一列
    energy = np.where(energy == 0,np.finfo(float).eps,energy) # if energy is zero, we get problems with log
    remove_ixs = np.where(energy <= 0.0001)[0]
    np.delete(pspec,remove_ixs,0)

    # np.where(condition, x, y)
    # 满足条件(condition)，输出x，不满足输出y。都有输出
    fb = get_filterbanks(nfilt,nfft,samplerate,lowfreq,highfreq)
    #np.dot()矩阵乘法
    feat = np.dot(pspec,fb.T) # compute the filterbank energies计算滤波器组能量，，，，，fb.T将数组转置
    feat = np.where(feat == 0,np.finfo(float).eps,feat) # if feat is zero, we get problems with log

    return feat,energy
def get_filterbanks(nfilt=20,nfft=512,samplerate=16000,lowfreq=0,highfreq=None):
    """Compute a Mel-filterbank. The filters are stored in the rows, the columns correspond
    to fft bins. The filters are returned as an array of size nfilt * (nfft/2 + 1)
    计算一个梅尔滤波器组。 筛选器存储在行中，列对应于fft bin。 过滤器以大小为nfilt *（nfft / 2 + 1）的数组形式返回

    :param nfilt: the number of filters in the filterbank, default 20.
    :param nfft: the FFT size. Default is 512.
    :param samplerate: the samplerate of the signal we are working with. Affects mel spacing.
    :param lowfreq: lowest band edge of mel filters, default 0 Hz
    :param highfreq: highest band edge of mel filters, default samplerate/2
    :returns: A np array of size nfilt * (nfft/2 + 1) containing filterbank. Each row holds 1 filter.
    一个大小为nfilt *（nfft / 2 + 1）的np数组，其中包含filterbank。 每行包含1个过滤器。
    """
    highfreq= highfreq or samplerate/2
    assert highfreq <= samplerate/2, "highfreq is greater than samplerate/2"

    # compute points evenly spaced in mels
    lowmel = hz2mel(lowfreq)
    highmel = hz2mel(highfreq)
    melpoints = np.linspace(lowmel,highmel,nfilt+2)
    # our points are in Hz, but we use fft bins, so we have to convert
    #  from Hz to fft bin number
    bin = np.floor((nfft+1)*mel2hz(melpoints)/samplerate)

    fbank = np.zeros([nfilt,nfft//2+1])
    for j in range(0,nfilt):
        for i in range(int(bin[j]), int(bin[j+1])):
            fbank[j,i] = (i - bin[j]) / (bin[j+1]-bin[j])
        for i in range(int(bin[j+1]), int(bin[j+2])):
            fbank[j,i] = (bin[j+2]-i) / (bin[j+2]-bin[j+1])
    return fbank

def hz2mel(hz):
    """Convert a value in Hertz to Mels

    :param hz: a value in Hz. This can also be a numpy array, conversion proceeds element-wise.
    :returns: a value in Mels. If an array was passed in, an identical sized array is returned.
    """
    return 2595 * np.log10(1+hz/700.)

def mel2hz(mel):
    """Convert a value in Mels to Hertz

    :param mel: a value in Mels. This can also be a numpy array, conversion proceeds element-wise.
    :returns: a value in Hertz. If an array was passed in, an identical sized array is returned.
    """
    return 700*(10**(mel/2595.0)-1)

def lifter(cepstra, L=22):
    """用一个倒谱提升器来提升倒谱矩阵。 这有增加高频 dct 系数幅值的作用.

    :param cepstra:  mel-cepstra 的矩阵，大小是 numframes * numcep。.
    :param L: the liftering coefficient to use. Default is 22. L <= 0 disables lifter.
               升降系数可以使用，默认值是22.10，升降器无法升降。
    """
    if L > 0:
        nframes,ncoeff = np.shape(cepstra)
        n = np.arange(ncoeff)
        lift = 1 + (L/2.)*np.sin(np.pi*n/L)
        return lift*cepstra
    else:
        # values of L <= 0, do nothing
        return cepstra

def delta(feat, N):
    """Compute delta features from a feature vector sequence.

    :param feat: A numpy array of size (NUMFRAMES by number of features) containing features. Each row holds 1 feature vector.
    :param N: For each frame, calculate delta features based on preceding and following N frames
    :returns: A numpy array of size (NUMFRAMES by number of features) containing delta features. Each row holds 1 delta feature vector.
    """
    if N < 1:
        raise ValueError('N must be an integer >= 1')
    NUMFRAMES = len(feat)
    denominator = 2 * sum([i**2 for i in range(1, N+1)])
    delta_feat = np.empty_like(feat)
    padded = np.pad(feat, ((N, N), (0, 0)), mode='edge')   # padded version of feat
    for t in range(NUMFRAMES):
        delta_feat[t] = np.dot(np.arange(-N, N+1), padded[t : t+2*N+1]) / denominator   # [t : t+2*N+1] == [(N+t)-N : (N+t)+N+1]
    return delta_feat

def read(data_path):
    '''读取语音信号
    '''
    # feat=[]
    # appendEnergy = True
    # if appendEnergy:
    #     feat[:, 0] = numpy.log(energy)
    wavepath = data_path
    f = wave.open(wavepath, 'rb')
    params = f.getparams()
    nchannels, sampwidth, framerate, nframes = params[:4]  # 声道数、量化位数、采样频率、采样点数
    str_data = f.readframes(nframes)  # 读取音频，字符串格式
    f.close()
    wavedata = np.fromstring(str_data, dtype=np.short)  # 将字符串转化为浮点型数据
    wavedata = wavedata * 1.0 / (max(abs(wavedata)))  # wave幅值归一化
    return wavedata, nframes, framerate
if __name__ == '__main__':
    data_path = r'H:\sharedir\file1\persion1\f0001_us_f0001_00001.wav'
    wavedata, nframes, framerate = read(data_path)
    feat = get_feature(wavedata)
    print(len(feat),len(feat[0]))#467 13